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Abstract. The leapfrog integrator is widely used because of its excellent 
stability in molecular dynamics simulation. This is recognized as being due to 
the existence of a discrete variational structure of the equations. We introduce a 
modified leapfrog method which includes an additional energy-like conservation 
law by embedding a molecular dynamics simulation within a larger dynamical 
system. 



1. Introduction 

The leapfrog integrator for molecular dynamics is known to display several exceptional 
features which allow it to have superior long-time dynamic stability compared 
with many higher-order integrators. This exceptional nature is due to the exact 
conservation of momentum and symplectic structure by the discretization (Hairer 
et al. 2002). The symplectic nature is particularly important for statistical mechanics 
applications because it tells us that there is a well defined density in phase space 
which is conserved as in the classic Liouville theorem; this density is the basis for 
the construction of the Gibbsian approach to statistical mechanics (Gibbs 2010). In 
addition the existence of a backward error analysis (Reich 1999) tells us that we are 
observing trajectories which are those of a perturbed Hamiltonian which is close to 
that which we are interested in. 

Despite these numerous advantages there is one important quantity which is not 
exactly conserved, this is the energy. It fiuctuates on small time scales and can drift in 
the very longest simulations. This is normally countered by introducing a coupling to a 
thermostat (Nose 1984, Hoover 1985, Frenkel & Smit 2002), but this leads to a change 
of ensemble from micro-canonical to canonical. The aim of the present paper is the 
construction of an integrator that is similar in many ways to the leapfrog method, but 
which is embedded in a larger dynamic system. The dynamics of the larger system 
are such that the original energy emerges as an additional conservation law. 

There are many trivial (and bad) manners to impose such an energy conservation- 
for instance one can regularly rescale the particle velocities. However, such arbitrary 
modifications to the dynamics break the symplectic structure which is a disaster for 
applications in statistical mechanics. Our approach to build a larger dynamic set of 
equations via variational methods in such a way that we can construct the discrete 
Hamiltonian of the system and thus explicitly understand the phase-space structure 
of the extended dynamical system. 
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This extended dynamical system is built using several components. We start by 
considering a discretized Lagrangian, in which the extra conservation is imposed by a 
Lagrange multiplier. From this dynamic system we build a discretized Hamiltonian. 
This Hamiltonian has a common defect that occurs in constrained systems (such as 
electrodynamics) - there is no momentum which corresponds to the multiplier. The 
solution is to add additional terms to the Hamiltonian which are zero, but which 
nevertheless generate independent dynamical equations for the multiplier. The logic 
is very close to the treatment of the potential in electrodynamics which has the formal 
role of being the Lagrange multiplier in imposing Gauss' law (Dirac 2001). 



2. Variational integrators 

We firstly resume how to pass from a discretized Lagrangian to the leapfrog integrator 
before generalizing to our more complicated constrained system: Newtons equations 
of motion for particles moving under velocity independent forces can be found by 
considering the variational problem 

^ I {^[^{ (1) 




In the following we will take all masses to be identical, and will allow q to denote a. Nxd 
dimensional vector corresponding to N particles moving in d dimensional space. This 
Lagrangian can be discretized by replacing derivatives by finite differences evaluated 
every r so that tk = kr: 

Lk = m — V{qk) (2) 

The discretized action principle is then (Guo et al. 2002, Marsden & Raiu 1999) 

5Y,rLk--Q (3) 

fc 

This variation then gives simple partial derivatives with respect to Xk so that 

m(gfe+i+gfe_i-2gfc) + y'(gffe) = (4) 

which is indeed a version of the leapfrog algorithm (Frenkel & Smit 2002). In the 
continuous time limit the energy is exactly conserved. 

However any time stepping procedure such as eq. ([4]) leads to a breakdown in the 
conservation of energy. The first step of our modified procedure will be to take the 
energy eq. (O and add it as a constraint to the original Lagrangian density eq. 
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Lk = m — V{q, 
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Afc(m^^^±l^ + yfe)-t/) (6) 

We will call the second line of eq. ([5]) the quasi-energy. In the following the dynamic 
schemes that we propose will conserve this discretized quasi-energy to machine 
accuracy. In this expression Afe is a Lagrange multiplier whose dynamics will be 
developed in the following sections. In the continuous time limit clearly A = 0. 
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The main questions that wih arrise are the following: In the presence of the 
extended dynamical system eq. © how do we interpret Liouville's theorem? What 
are the corresponding momentum variables for the discretized evolution equations that 
comes from eq. ([B]). In order to answer these questions we pass from the Lagrangian 
description to a Hamiltonian form for the dynamics. 



3. Discrete Hamiltonians 

We will need to introduce a slightly more formal notation in order to pass from the 
above Lagrangian formulation to a discrete Hamiltonian form. However this notation 
is such that we find expressions which are very close to those in standard treatments 
of Hamiltonian dynamics. We firstly introduce the finite time difference operator 

Agfc = {qk+i -qk)/T (7) 

We also need its adjoint A* which is defined so that 

^flfeAgfe = - '^qk{ak - ak-ij/r = ^afeA*^^ (8) 

Thus we see that 

A*gfc = ^^qk-i (9) 

There is a natural shift of unity in indices when performing the discrete analogy of 
integrating by parts. 

The Lagrangian equations of motion are then 

A*^ + |^^^0 (10) 
9Agfe dqk 

which is very close to their form in the continuum limit. We now define the momentum 
variables: 

dL 

^'fe+i = flA — = m((7fe+i - gffc)(l + Afc)/r (11) 
and construct the Hamiltonian as usual 

H{qk,Pk+i) = Pk+i^qk - L 

= .(/h\ +Vk{l-Xk) + XkU (12) 
2(1 + AfejTO 

Let us consider the equations of motion which come from applying Hamilton's 
principle to eq. ([T2|) . We calculate 



Sy\pk+iAqk-Hipk+i,qk)]^0 (13) 



k 



and find 



^ = -Apk (14) 
oqk 

dH , 

7^ = Agfe (15) 

OPk+l 

which is the discretized form of the Hamiltonian equations of motion. More explicitly 
we have 

Apk = Pk+i - Pk = -Til - XkM (16) 

Agfe = qk+i -qk = t ^"^^^ (17) 
m(l + Afe) 
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When Afc = this corresponds to the standard ahernating update in the leapfrog 
algorithm. 

Now consider the equation which comes from varying the Lagrange multiplier : 

l^^^-^('^-)- 2..a + A..)- =^--» 

which is just the equation for the quasi-energy in the Hamiltonian picture. Note 
this quasi-energy conservation does not imply that the Hamiltonian eq. (jl2p is itself 
conserved, however if they are numerically close one might hope that the stability of 
the algorithm is also improved for H. 

If Afc were a true dynamic degree of freedom we would have deduced from eq. (jl8p . 
in analogy to eq. dH]) that 

-A7rk = Wk{qk,Pk+i,Xk) (19) 
where tt^ the conjugate momentum to Afc. We discover that in order to have a full 
Hamiltonian description of the system we must add this extra degree of freedom, but 
also that tt^ = for all k in order conserve the quasi-energy. We will show this is 
possible later, but firstly move on to the practical question of implementation of the 
algorithm. 

4. Integration loop 

The equations eq. eq. (|17p together with the constraint equation Wk = 0, eq. p^ . 
tell us how the positions and momenta of the particles evolve within a time step. We 
now show that the equations have explicit (non-iterative) solutions that require only 
small modifications of the usual leapfrog step. 

We firstly take eq. ([TSl) and express Pk+i in terms of pk- 



o nr ^/^ (Pk + Tjl - \k)fkY 

2'"^^ - = WTXkf ^^^^ 

or 

S{l + \kf ^pl + 2Tpkfk{l~Xk)+T''fl{l-\kf (21) 
with fk the force. Thus 

\1{S - T^fl) + \k{2S + 2Tpk ■ fk + 2t2/|) + 

{S-2Tpk- fk-r\fl-pl) = Q (22) 

with S = 2m{U — Vk). This is a simple quadratic equation for A^ which involves 
quantities which are already calculated within a leapfrog integration loop. In practice 
A remains small throughout our simulations and its value is close to 

A « -rp • f/S (23) 

We can integrate the equations with the following loop 

know Xk-i,qk,Pk 

calculate fkiqk) 

calculate Afc eq. ((22|) 

calculate Pk+i eq. ITB)) 

eq. (fT5)) satisfied at this moment of cycle 

calculate qk+i eq. ([T71) 

know Xk,qk+i,Pk+i 
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This is the practical generahsation of the generahsed leapfrog method with the 
addition of an exact conservation of the quasi-energy. 

5. Momentum for A 

The equations of motion as stated above are adequate to implement the algorithm. 
However they contain a formal weakness. While p and q evolution is the result of 
a variational principle eq. (|13|) this is not true of A. The Lagrange multiplier is 
simply slaved to impose energy conservation. Such slaved variables are well known 
in quantum mechanics, indeed the electrostatic potential is an example of such a 
variable; this explains the initial difficulties in the quantisation in electrodynamics 
due to the lack of an obvious conjugate momentum. We now show how to render the 
equation for the evolution of A autonomous, and thus better understand the phase 
space of the enlarged dynamic system. We will use methods which are rather similar 
to those invented in electrodynamics (Dirac 2001, Leimkuhler & Skeel 1994) where 
the electrostatic potential also has a role which is similar to a Lagrange multiplier. 
There are clear analogies too with our previous work on local simulation algorithms 
for charged media (Rottler & Maggs 2004, Maggs 2002). The first problem with 
the equation eq. ()12|) is that it does not include a momentum variable, iTk which is 
conjugate to A^. We correct this deficit with the following ansatz: we add an extra 
term 

7rfc+iAifc((jfe,pfc+i, Afc) (24) 

where fik is a function that we will construct later. This leads to the following 
equations of motion: 

OH 

AAfc = = t^k (25) 

dH dill. , , 

oAfc (9Afc 
We now use the idea of weak constraints: For arbitrary functions /i^ we 
have a Hamiltonian system. We will show that the correct choice of the function 
Mfc(pfe+ij 9fci '^fe) allows us to impose both the conservation of eq. (|18|) but is also 
compatible with Avrfe = 0. We then start the dynamic system in the state ttq = and 
this remains true for all further times in the dynamics. 

We now show that the function /x^ is indeed only a function of the objects 
Afc). We proceed by considering W-k at two successive time steps, imposing 
conservation of the quasi-energy 

W^fc((7fc,pfc+i, Afc) = tyfc+i(gfe+i,pfc+2, Afc+i) (27) 

We now eliminate the variables (7fc+i and pfc+2 from the right hand side using eq. ()17p 
and eq. ([T6l) . If we do so the right hand side of eq. (|271) is a function of g/cp/c+i, Afc+i 
We can thus, in principle solve for A^+i and write the evolution in the form 

Afc+i = Afc + fikiqk,Pk+iAk) (28) 

as needed for the dynamics of iTk , eq. ([M)) , eq. (|25l) 
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6. Phase space and Liouville 

We have succeeded in embedding our original system of particle dynamics in a larger 
system such that the quasi-energy is exactly conserved. To do so we were obliged to 
introduce two new variables Afc which started as a simple Lagrange multiplier and tt^ 
which is the conjugate momentum. For a system of N particles in d-dimcnsional space 
this gives us a phase-space of dimensions 2dN + 2. 

We now study the Jacobian of the discrete evolution equations to show that 
the constrained dynamical systems are compatible with the assumed measure. The 
maps eq. and eq. (fT7)) are easily seen to have unit Jacobians on the phase space 
defined by the variables {q,p, A, tt). The evolution equations for A require slightly more 
study. The Jacobian for A^+i = A/c + Hkiqk,Pk+iAk) is given by J = 1 + d^k/dXk- 
while we see we can re-arrange eq. ((26|) to give TTk+i = 7rA_,/(l -I- d^k/dXk)- It is thus 
the product of these two factors which ensures that the Jacobian of a full time-step 
is indeed unity. We thus see that introduction of the "dummy" momentum tt has 
absorbed the fluctuations in phase space volumes that would otherwise result from 
the use of the Lagrange multiplier. 

Thus we have a complete set of Hamiltonian dynamics on the extended phase 
space with the extra pair of variables A, tt which are now fully autonomous. We 
however choose special initial conditions tt = that lead to the exact imposition 
of the energy conservation. We conclude that we have a phase space measure of 
the form dq dp dX dn, with conservation laws imposing constant quasi-energy, particle 
momentum, and n (Khinchin 1949). 

7. Time reversal 

While we have added a conservation law to the leapfrog integrator we have also lost a 
symmetry which is present in the standard leapfrog algorithm - it is time reversible: 
The more complicated quasi-energy conserving version does not re-trace its trajectory 
when momenta are reversed. This extra symmetry can be imposed in our algorithm 
by alternating the direct step (described above in section 4) with a version in which 
each step in implemented in reverse order, with t ^ —t. The main technical difficulty 
is that the equation of Xk-i given qk and pk becomes implicit and must be solved by 
iteration. In practice we find that a simple iteration procedure converges to machine 
precision in two steps if we use eq. (|23p as a starting guess for A. The algorithm 
in which a the direct and reversed step are alternated then displays time reversal 
symmetry. 

8. Conclusion 

We have constructed a variational integrator which includes an additional conserved 
quasi-energy. Due to the variational Hamiltonian form we are able to study the 
phase space measure and understand the discrete Liouville theorem that is implied 
by the dynamics. Implementation of the algorithm requires a small overhead 
in computational effort compared with the standard leapfrog integrator. We 
have implemented a version of the code for the molecular dynamics study of a 
truncated Lennard- Jones potential and verified the stability of the quasi-energy during 
simulation. 
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